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Abstract 

Sampling from flat energy or density distributions has proven useful in equilibrating 
complex systems with large energy barriers. Several thermostats and barostats are 
presented to sample these flat distributions by molecular dynamics. These methods use a 
variable temperature or pressure that is updated on the fly in the thermodynamic 
controller. These methods are illustrated on a Lennard- Jones system and a structure- 
based model of proteins. 



I. Introduction 

Molecular dynamics (MD) is a useful tool to study complex molecular systems. 
A regular MD simply follows Newton's equation, which conserves the total energy E and 
thus generates configurations in a microcanonical ensemble. We can, however, modify 
the equation of motion by an artificial thermostat 1 " 9 to sample a different distribution, in 
which the total energy fluctuates, such as the canonical distribution. 

For a system with high energy barriers, the limited exploration of the energy 
landscape in the canonical distribution is insufficient. A multicanonical ensemble 10 " 13 , in 
which the energy is broadly sampled, is more suitable, see Fig. 1(a). The construction of 
the multicanonical ensemble is nontrivial, for the flat energy distribution requires the 
inverse of the unknown density of states, \/D.(E) , as the sampling weight. 

We here pursue the construction of multicanonical 10 " 13 MD by using a variable 
temperature equal to f3{E) = dlogD./dE . We describe several thermostats to set this 
temperature. The method is built on a previous Monte Carlo (MC) method by Yan and 
de Pablo 14 . Here the sampling along the total energy is directly incorporated into the 
variable-temperature thermostat for the MD simulation. 

In section II, we describe the method. In section III, we numerically test the 
method in a few examples. We conclude the article in section IV. 



II. Method 



II.A Microcanonical, canonical, and multicanonical ensembles 
Microcanonical ensemble 

In classical statistical mechanics, integration of Newton's equation generates 
configurations in a microcanonical ensemble under a fixed total energy 
E = U(q) + K(p) , where U (q) and K(p) = p 2 /2m are the potential and kinetic energies, 
respectively. Thus, configurations from an ergodic trajectory sum to the density of states: 



where 0(x) is the step function, which is 1 if x > or otherwise, and S(x) = &(x) is 
the ^-function. We have also integrated out the N f momentum degrees of freedom 15 : 



with T(x) being the gamma function. For a molecular system, N f is three times the 

number of particles less the number of conserved quantities, such as the total linear and 
angular momenta. 

The microcanonical temperature, or the derivative of the entropy 
S{E) = k B logQ(£) 16 ' 17 with k B being the Boltzmann factor, can be evaluated as a 
configuration average 15 ' 18 



Q(£)ocf S[E-U(q)-K(p)]dqdpoz f [E -U(q)f f/2 l ®[E -U(q)]dq, (1) 
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where f3{E) = \/[k B T(E)] is the inverse temperature, (...) £ denotes an average under a 



fixed total energy E, and we have used Eq. (1): 



1 dQ(E) = r ( N,/2 -l ) [E- U( q )] N '/ 21 ©[E - U(q)] = ( N,/2-l \ 
Q(E) dE **{E-U(q)) Q(£) \ K J E 

Here we have assumed N f > 2 so that the singularity at E = U (q) is negligible. Other 

18 

expressions for /3{E) also exist , as shown in Appendix A. 



Canonical ensemble 

The canonical ensemble is a superposition of microcanonical ensembles at various 
E , each with a weight exp(-/3E) . The overall E distribution is given by 

p(E) oc f S[E -//(q,p)]Jq dp exp(-/?£)= Q(£)exp(-/?£) . (3) 

To realize a canonical distribution in MD, one needs a thermostat 1 " 9 to alter the equation 
of motion either deterministically " or stochastically 



Multicanonical ensemble 

We now consider a generalized ensemble composed of microcanonical ensembles 
with an arbitrary weight exp[-/ (E)] instead of the Boltzmann weight exp(-/3E) . We 

approximate / (E) by a piecewise linear function on a finely-spaced grid E ,E v ...,E n , 

as in Fig. 1(b). Then within each (E k ,E k+l ), f(E) « f(E k ) + f'(E k )(E -E k ) , and the 

ensemble can be treated as locally canonical with the effective temperature being 
f'(E k ) 10 . The constant f(E k ) - f'(E k )E k does not distinguish configurations, and is 

therefore irrelevant. In the limit of infinitesimal intervals, the local effective temperature 
becomes f'(E) . 



The multicanonical ensemble is the special case of the generalized ensemble 
described above, in which / (E) = logQ(is) + const. , such that the E distribution 

An»c,(£)= f S[E - //(q,p)] dq dp exp[-/(£)]oc Q(£)exp[-logQ(£)] = 1 
J <\,p 

is flat. Thus, the local temperature f'(E) should be energy dependent, and it is equal to 
(logQ)'(E) , which can be evaluated by Eq. (2). 

Thus, molecular dynamics in the multicanonical ensemble can be achieved by a 
thermostat with an energy-dependent target temperature f'(E) = (logQ)'(is) . Below, we 

describe how this can be done. As in the previous MC method 14 , Q(E) is estimated from 
its derivative by Eq. (2), instead of from the energy histogram 10 ' 12 ' 13 . Strictly speaking, 
detailed balance is violated if we constantly update (logQ)'(E) by Eq. (2) in simulation. 

Earl and Deem, however, showed that this practice, commonly used in the adaptive 
integration methods 14 ' 19 and to be used here, satisfies balance in infinite-order Markov 
chains 20 . This updating practice converges in the present application, and so the present 
method satisfies balance at long times. Note also that the multicanonical ensemble here 
is formulated for the total energy 14 , instead of the potential energy 10 ' 12 . 

II.B Algorithm 

We now describe an algorithm for a multicanonical MD. The aim is to sample a 
flat energy distribution in a given range (E mill ,E max ) ■ We split the range into n intervals 

at E = E mia , E x , E n = E max , such that each (E k , E k+X ) accumulates data for 

computing the average of (N f /2 - 1)/K(p) to estimate /?(£) in Eq. (2): 



The algorithm: 

1 . Integrate Newton's equation for At MD . 

2. Locate the interval k that contains the current total energy E (if E min < E < E max ); 
and deposit (N f /2-1) / K(p) into the accumulator there. 

3. Use Eq. (2') to estimate fi(E) ; but if the number of data points is less than 100, or 
if E lies outside of (E min , E maK ), use a default value /3(E) = /3 & , which should 
correspond to an energy E e (E min ,E mRX ) . 

4. Use fi{E) in one of the variable-temperature thermostats listed below. 

5. Repeat Step 1 until the simulation ends. 

Remarks 

1 . We show in Appendix A that the use of the interval-averaged formula Eq. (2') to 
approximate /3{E) produces no systematic error: the estimated values of the 

density of states \ogQ.{E k )= I /3{E)dE at the boundaries E k asymptotically 

approach \ogQ,{E k ) for any interval size E k+l -E k . 

2. While the average in Eq. (2') is invariant with respect to the interval size, the 
energy range, the number of intervals, and the number of data points will affect 
the rate of convergence. The number of intervals should be large enough so that 

the logQ(is) is sufficiently close to logQ(is) , i.e., the dotted line segments in 



Fig. 1(b) closely approximate the solid curve. Furthermore, the number of data 
points per interval should be sufficiently large to provide a reliable estimate of the 
average. 

3. In Step 3, a constant fi(E) = J3 A is assumed for E outside of (E min ,E max ) for 

simplicity. It, however, causes a large discontinuity there; and a better boundary 
condition may improve this. The number 100 is also arbitrary. 

4. After simulation, we can recover the density of potential-energy state 

g(U) = \ S[U - U (q)]dq from the observed distribution h(U) of the potential 

Jq 

energy Uby 



g(U) = h(U)l 



\™{E-U) Nf l 2 - l l£l{E)dE 

see Appendix B for a proof. The denominator of Eq. (4) can be evaluated exactly 
by the incomplete gamma function 21 , as logQ(ii) is piecewise linear. One can 
then compute the canonical partition function Z{/3) = j g(U) exp(-/? U) dU , and 

its derivatives. Eq. (4) also applies to a Q(E) that does not produce a flat E 
distribution, see examples in Appendix A. 



(4) 



II.C Variable-Temperature thermostats 

The thermostats for a canonical ensemble need to be modified as follows for a 
multicanonical MD. We list a few examples below. The dimension is denoted by d , 
which is usually 3. 



Nose-Hoover chain thermostat 

The Nose-Hoover chain equations 5 ' for a canonical ensemble are 

q = p/m, p=F-C lP , ^={v 2 lni-N f l/3)lQ x -^, 

t k = {QaJ-VP)lQ k -^U (for£ = 2,...,M-l) (5) 

L={Qm-£mJ-vp)IQm, 

where £ k are the thermostat variables, and Q k their masses. The corresponding 
Liouville's equation 

dp/dt + d(qp)/dq + a(pp)/e P + Ylji^P)l^k = o (6) 
has a stationary (dp/dt = 0) solution 

p(q,p,{CJ) oc exp[-/?£/(q) - /?p 2 /2m - fi^LO* til 7 * ' 
which is indeed canonical 5 . 

For a multicanonical ensemble, we change the equation for £ in Eqs. (5) to 

£ = {[0(E)/ fi ] \> 2 /m - N f lP }/Q x - , (5') 

and we replace the /? in the equations for £ k , k > 2 , by J3 , such that the new stationary 
solution of Eq. (6) is 

p(q, p, } ) oc exp {- log Q [U(q) + p 2 /2m] - & ^ Q k fa , (7) 
where /? is a reference inverse temperature, against which the set point of the system's 
inverse temperature is scaled in the term in braces in Eq. (5'). That is, with the use of 
/3(E) in Eq. (5'), the stationary distribution is the desired Eq. (7), no matter the value of 
/? . One expects that /3 Q values near the average of /3(E) will be efficient for sampling. 



We recover the original Nose-Hoover thermostat by setting £" 2 = , and removing 
the equations for ^ with k > 2 and for C, M . 

Velocity-Rescaling thermostat 

The velocity-rescaling thermostat 9 achieves a canonical distribution by re- 
sampling the distribution of the total energy E under fixed q coordinates: 

p(E)dE = [E- U(q)f f/2 ~ l Qxp(-fiE)dE , 
cf. Eqs. (1) and (3). The change E — > E' , hence K — > K' , is realized by scaling the 
velocity or momentum: p — > Jk' I K\) , for U (q) is fixed. The kinetic energy is 
randomly updated and follows the Langevin equation 9 

dK/dt = \N f /p -K + pK/j3 £ , (8) 

where £, is a Gaussian white noise that satisfies (£(0£(O) = 8(t-t') . 
For a multicanonical ensemble, we modify Eq. (8) to 

dK/dt = \N f /fi -iME)//3 ]K + pK/fi 4 , (8') 
such that the corresponding Fokker-Planck equation 

dp/dt + d({±N f -[ME)//3 ]K}p)/dK-d 2 (Kp//3 )/8K 2 =0 

has the desired stationary solution p(K) cc K Nf ' 2 l /Q.(U + K) . Here, f3 a is a reference 
inverse temperature. 

Note that Eq. (8) or (8') can be integrated with a different time step from the MD 

one. 



Monte-Carlo velocity-rescaling thermostat 

We can also replace Eq. (8') by an explicit MC move in the E, or K space 14 . We 
first propose a change of log^ by 8 e (-A, A) , then accept the new K' = Ke s by 



A(K -> K ) = mm\ 1. 



K 



Q(K' + U) 

If iT' is accepted, it is also realized as a velocity rescaling: p — > ^K' / Kp . 

Langevin dynamics 

Langevin dynamics 1 ' 17 modifies Newton's equation of motion by introducing a 
velocity damping term and a thermal noise: 

q = p/m, i = F-C V + j2m£/fi§, (9) 

where C, is a damping constant, and £ is a Gaussian white noise satisfying 

(zm j (t , ))=$(t-t')s ij . 

For a multicanonical ensemble, we change Eq. (9) to 

q = p/m, p = F - [A£)//WP + pmZ\p& (9') 
such that the corresponding Fokker-Planck equation 

dp/dt + d[( V /m)p]/dq + 8{[¥- P(E)//] £p]p}/dp - d 2 [(mC / /3 )p]/dp 2 = 

has the desired stationary (dp/dt = 0) solution p(q,p) <x l/Q[{7(q) + p 2 /2m]. Here, /3 
is a reference inverse temperature. 



Andersen thermostat 

In the Andersen thermostat 2 ' 8 , we replace the velocity or momentum p° ld of a 
random particle i by a Gaussian random vector drawn from the Maxwell distribution: 

/Wen(# PD = (P/2xmf 2 Q xp[-/3(prf/2m] . (10) 

For a multicanonical ensemble, we treat p° ew from Eq. (10), with /? replaced by 
/3(E) , as a MC trial, and accept it by the Metropolis probability 

^(pf d ^ P r) = minjl, ^ {E) ^M a xweii(A£');P°' d ) l (n) 

where E' = E + (p; ew ) 2 /2m - (pf d ) 2 /2m , and the two Maxwell factors /? Maxwell (. . .) , 
defined in Eq. (10), offset the a priori sampling bias. For the canonical ensemble, 
(3(E) = (3(E') = (3 , Q(E)/Q(E') = Qxp[/3(pf d ) 2 /2m-/3(tf ew ) 2 /2m], and Eq. (11) yields 
unity so that the trial is always accepted. 

II.D Volume distribution 

We can extend the above sampling strategy to explore a broad volume 

22 23 

distribution, which is useful in studying a liquid-gas phase transition ' . A flat volume 
distribution can be inefficient, however, as drastic free energy changes occur only at a 
high-density, or small-volume, regime. We therefore sample the distribution 
p(V) oc l/V a . To sample a flat density distribution, p n {n) = const, with n = N/V , we 
use a = 2 , because 

p{V) dV = p n (n) \dn\ = p n {n)N dV/V 2 oc (\/v 2 ) dV . 



To sample the distribution p(V)x V a , we replace the constant target pressure in 

2 5 7 8 24 26 

an isothermal-isobaric simulation • • > > ~ by the variable 

P(V) = ~ 7— r (12) 

V ) 

\ /^-<K<K t+1 

where /^.foP) = (p7 w + F • q)/(<* ■ V)-d(/) m {V)ldV , where $ ail (7) is the energy tail 

7 8 24 26 

correction, such that the potential energy U (q,V) = U tmriC (q) + <fi t!lil (V) ' ' " . The weight 



V" in Eq. (12) ensures unbiased free energy F(V k ) = -I p(V)dV for a finite interval 
width V k+l - V k , see Appendix A. 



Nose-Hoover chain barostat 

To sample p(V)cc V~ a , we modify the Nose-Hoover chain 5 ' 7 ' 8 equation for an 

isothermal-isobaric ensemble to 

q = p/m + ?7q, p=F-(^ + ^)p, V = dr/V, 
V = {[A„,(q,P) " P(VW + (1 - a)lfi}d/W - £ 7 , 

£ = [p 2 /™ + JF/; 2 - (tf, + l)//?]/Q - , (13) 
£ = (QUJ ~ Vfl/Ot -CwCt, (for * = 2, . . . ,M - 1) 

c m =(q m - 1 CmJ-v^/Qm, 

where ^ and ^ are the variables for the barostat and thermostat, respectively, with W 
and Q k being their respective masses. The corresponding Liouville's equation 

dp/dt + d(qp)/dq + d(j>p)/dp + d{i)p)ld V + Z M k J(C k p)/dC k + d(Vp)/8V = 



has the desired stationary solution: 

p(q,p,7,{CJ,F)^exp{- y 9[t/(q,F)+pV2m + r^ 2 /2 + ^ i a^ 2 /2-^)]}^ 



rV 

where F(V) = -I p{V')dV is the estimated Helmholtz free energy. The term I- a in 

the equation for 77 accounts for the ensemble weight V" . We recover an isothermal- 
isobaric sampling by setting p{V) to a constant and a — > . The 1 - « — > 1 corrects the 
small error in the original Nose-Hoover barostat 5 ' 7 ' 8 . 

We recover the Nose-Hoover case by setting £" 2 = , and removing the equations 

for £ k with k > 2 and for C, u . 
Langevin barostat 

Alternatively, the pressure can be controlled by a Langevin equation separate 
from the thermostat. The Langevin barostat is a stochastic version of Eq. (13): 

q = p/m + /7q, p=F-^p, V = dr]V, 

V = {[A„,(q,P) " P(VW + (1 - a)lf3}d/W + pC/(0W)Z. 
where S, is a Gaussian white noise satisfying (i;(t)i;(t')) = S(t-t') , g and W are two 
constants. The Fokker-Planck equation 

dp/dt + d(qp)/8q + d(pp)/8p + 8(Vp)/dV + d[(r/) det p]/dr] - d\{^ 'I J3W) p]/ 'drf = , 
where (r/) det = {[p int .(H,p) ~ p(V)W + (1 _ a)/ P)d/W - £r], admits the desired stationary 
distribution: 

p(q,p,r ? X)^^p{-fi[U(qX)+pV2m+Wr ? 2 /2-F(V)]}V- a . 
Explicit Monte-Carlo sampling 

8 22 

The above Langevin equation can be replaced by a MC move ' : we first propose 
a change of \ogV by 8 e (-A, A) , then accept the new V = Ve by 



A{V V) = min\ 1, 



exp {-p [Uijjr/Vq, V) + (^/^p) 2 1 2m - F{V')] } V 
exp {-£ [t/(q, + p 2 /2m - F{V)] } V~ a 




. (15) 



If the move is accepted, we change the coordinates as q — > yV'/V q , and momenta as 
p — > ijV/V'p . The phase-space volume element dqdp is conserved. 

III. Examples 

III.A Flat energy distribution 

We tested the method in the jV = 108 particle Lennard- Jones system 8 . The 
potential was cut off at r c = 2.5 A and shifted to zero. The density was n = N/V = 0.7 . 

Each trajectory was integrated for 5 x 10 8 steps; and the time step was At MD = 0.001 
unless specified otherwise. The energy E range was {-AN, 2N) , and the interval was 
0.027V . The default inverse temperature f3 A for E outside of the range or insufficient 
data and reference inverse temperature /? were both 1.0. 

The energy distributions p{E) from simulations using several variable- 
temperature thermostats are shown in Figs. 2(a)-(e). The flatness of the energy 
distributions p{E) helped detect if the variable-temperature thermostat was correct. The 
modifications in Sec. II. C were necessary to produce flat distributions. That is, if the (5 

was simply replaced by f3{E) in Eqs. (5), (8), (9), and (10), flat energy distributions were 
not achieved, as shown in Fig. 2(a), 2(c), 2(e), and 2(f), respectively. 

The thermostat equations are stiff if the mass Q x in the Nose-Hoover chain is not 

extensive. An efficient Q x should be proportional to the number N f of degrees of 



freedom 4 ' 6 ' 25 , which can be seen from the scales of the equation for £ in Eq. (5'). From 
the left side, we have p(£ x ) °c exp(-/? Q x ^ x /l) by Eq. (7). Thus, £ oc \j ^Q x and 
£ ~ Ci/^md ~ V( a ^mdVS) • Fr o m th e ri g ht side of Eq. (5'), 

[P(E)I p ] p 2 /m-N f /P ~ ^iVy , since it is proportional to the standard deviation of 
p 2 /m . Thus, the right side ~ f j Q x . Equating the two results, we find 

/At MD ~ ^N~ f . The other Q k , k > 2 , should be of order unity 6 ' 25 . With g, =10, 
which was much smaller than , visible deviations from the flat distribution can be 
seen, but a smaller MD time step At MD recovered a flat distribution, see Fig. 2(f), as it 
must by Eq. (7). 

The approximate 15 temperature formula P(E) = N f j(2K} 14 ' 27 did not produce a 

flat energy distribution, but rather a p{E) cc , as shown in Fig. 2(d). This is 

because P(E) corresponds to a modified density of states weighted by 2K , as discussed 
in Appendix A. 

III.B Structure-Based protein model 

As an application of multicanonical MD, we studied the protein villin headpiece, 
Protein DataBank code: 1 VII 28 ' 29 , by a simplified structure -based model 30 " 32 . Details of 
the model are described in Appendix C. Each trajectory was run for 2 x 10 9 steps with 
A^ MD = 0.002 . We used the MC velocity-rescaling variable-temperature thermostat with 
A = 0.05 . We show below that the thermodynamics of the model depends critically on 



the cutoff distance r c that defines contacts; the interaction is attractive between contact 

atoms, but repulsive otherwise. 

As shown in Fig. 3(a), the microcanonical temperature functions /3(E) under 

several different r c were similar to pressure -volume isotherms around a liquid-gas 
transition in a simple fluid 17 ' 33 . The back-bending S-loop 34 that signifies a first-order 
phase transition occurred only when r c > 7 A. The sensitivity to r c may explain that both 

the presence and absence of S-loop were previously observed in similar models 13 ' 35 . A 
strong first-order transition in the energy space may be absent with the common 
definitions of contacts 31 ' 32 ' 36 , whose equivalent r c were about 4 to 5A. 

For r c > 6.5 A, the energy distributions p(E) around the transition temperature 
P* E were bimodal. Here, 0* E was defined as the temperature at which the two modes 
share the same height. See Fig. 3(c) for the case of r c = 8 A. The free energy profile 

- logQ(£) + /3* E E = f[/3* E - fi(E')]dE' indeed had two basins. 

We computed the density of potential-energy state g(U) by the reweighting 
formula Eq. (4). The transition temperature /?* determined from g(U) was similar to 
j3* E , see Fig. 3(d). The free energy profile - \ogg(U) + P* V U along U showed a higher 
barrier than that along - logQ(ii) + /?* E . 

III.C Flat density distribution 

We tested the extension in Sec. II. D for the flat density distribution also on the 
108 particle Lennard- Jones system. The pair potential was cut off at half-box length to 



facilitate volume change moves . The density range was (0.05, 0.75) with a spacing 
0.002. The default pressure p A for a density outside of the range or insufficient data was 

0.1. T = 12. Each trajectory was run for 4 x 10 8 steps with At MD = 0.001 . 

The density distribution, pressure-density curves, and free energy profile are 
shown in Fig. 4 for the Nose-Hoover chain Eq. (13), Langevin Eq. (14), and MC Eq. (15) 
barostats. The density distributions were flat. The obtained transition pressures were 
0.0696, 0.0695 and 0.0697 for the Nose-Hoover chain, Langevin and MC barostats, 
respectively. The free energy profiles obtained at the transition pressures also agreed 
with one another, see Fig. 4(b). 

IV. Discussion and Conclusion 

A multicanonical, or flat-energy-distribution, molecular dynamics can be 
conveniently realized by a variable -temperature thermostat. Similarly, a flat-density- 
distribution MD can be realized by a variable-pressure barostat. The method is easy to 
implement with minimal modifications from a regular constant temperature or volume 
simulation. We have demonstrated the method of a Lennard- Jones system and a 
structure-based model of the villin headpiece. 

The method computes the density of states Q(E) 14 instead of the density of 
potential-energy states g(U) 10 ' 12 . The reweighting formula Eq. (4) is used to recover 
g(U) . In this way, thermal averages in the canonical ensemble are computed after 
simulation. 

The method is somewhat more robust, in terms of achieving a flat energy 
distribution, with a smaller energy interval AE = E k+l - E k , although it is formally correct 



for any AE , as shown in Appendix A. A smaller interval also yields a finer Q.(E k ) , but 

may slow down the on the fly construction of the ensemble weight. The integral-identity 
technique 29 ' 37 , which allows borrowing data from neighboring intervals, and the adaptive 
averaging technique from the continuous tempering 29 may help improve convergence. 

The multicanonical ensemble is commonly used in Monte Carlo and molecular 
dynamics simulations to overcome energy barriers. We here have described a number of 
variable-temperature thermostats and variable -pressure barostats to achieve flat energy or 
density distributions. We expect these methods may prove useful in simulating systems 
with complex energy landscapes. 
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Appendix A. Unbiased density of states at interval boundaries 

We show here that the values of the density of states 0.(E k ) = I /3{E)dE at the 

boundaries points E k are unbiased, although fi(E) is an approximate average of /3(E) 
over finite intervals (E k ,E k+l ) . This is because 



Q(E k+1 ) Q(E k ) _ ce m d 

h(E k+l ) h(E k ) k dE 



n(E) 

Q(£) 



dE 



= r + ' [(logQ)'(£) - (logQ)'(£)]2^£ 
JE * Q,(E) 

'NJ2-X 



= K 



Me)) 



E t <E<E M 



where we have used (logQ)'(£) = [(N f /2-l)/K(p)) E , (logQ)'(E) = p(E) . The ratio 
Q(E)/q.(E) = I S[E - H(q,p)]dqdp Jo,{E) is the distribution density of the overall E 
distribution, and we have defined h k = f w Cl(E)/h(E)dE . Since the right hand side is 

JE k 

zero by Eq. (2'), the left side vanishes. So Cl(E k+1 )/Q.(E k ) = D.(E k+l )/D.(E k ) , and by 
recursion, Q(E k )/Q(E ) = Q(E k )/Q(E ). Thus, Q(E k ) differs from Q(E k ) only by a 



multiple. 



The alternative formulas for /3{E) : /? (p) (£) = N f /(2K(p)) ] and 



E„<E<E t+1 



in ref. 14,27, have a similar property, although 



they do not produce flat E distributions, even with infinitely small intervals. To see this, 
consider an arbitrary vector field v = v(p,q) . We define 

Cl D (E) = |(v • VH) S[E- H(q,p)]dq dp = Q(£)(v • VH 
Q N (E) = |(V • v) S[E - H(q,p)]dqdp = Q(E)(V ■ y) £ . 

Now integration by parts yields dd D (E)/dE = Q N (E) 1 ; and 



Cl D (E k+l ) Q D (E k ) _^ 
Q(E k+l ) Q(E k ) 



Q(£) 



(logQ)W 



=r 

JE 

= h ^-v) Ek<E<EM -kE)(WH) Ek 



Cl(E) 



Q(E) 



Q(E) 



dE 



<E<E t 



Thus, if we set fi(E) = (V • v) £<£<£ + _ /(v • VH)^^ , Q(E k ) will asymptotically 
approach the modified density of states Q. D (E k ) instead of 0,(E k ) . Accordingly, the 
converged energy distribution p{E) cc Q.(E)/Q. D (E) = \j(y ■ Vi/) £ is not necessarily flat, 
and the resulting Q,(E) should be divided by (v • Vi/) in order to recover Q,(E) . 

The two alternative mean force formulas mentioned above correspond to v (p) = p 
with v (p) -VH = v 2 /m, and v (F) = F with v (F) • WH = -F • F , respectively. Thus, the 
overall energy distributions p (v) (E) oc l/^p 2 /m^ and p (F) (E) oc l/(F ■ F) £ are not flat. 

We can force a flat distribution by using a "normalized" vector field. For an 
arbitrary v , we construct v' = v/(v • VH) , then 18 

yields the unbiased and a flat distribution, since (v' • Vi/) = 1 , and 

0, D {E) = D.(E) . The above two examples for v (p) = p and v (F) = F , after the 
normalization, become Eq. (2') and 3(E) = (V • [F/(F • VC/)]\ B B 18 , respectively. 

\ / E k <E<E t+l 



We can show Eq. (12) in a similar way. We define the canonical partition 
function at a fixed volume V as 



Q(V) = exp[-/JF(F)] = f exp[-/jtf F (q,p)]tfq dp , 

and its estimate 

log Q(V) = -pF(V) = p\ V p(V')dV' . 

Then, 



Q(V k+l ) Q(V k ) _ d 



= P + — W± dV = -p\ Vk+ \F\V)-F'{V)T 
Q(V k+1 ) Q(V k ) ^ dV[_Q(V)\ * 

= -fit' [FW) - F'(V)]V a w(V)dV 

= PK([ P (V)-p{vw a ) Vk<v 



Q{V) dV 

Q{V) V a 



<v<r. 



where w(V) = f exp[-/?// K (q,p)] dq dp /[Q(V)V a ] = Q(V)/[Q(V)V a ] is the overall 
volume distribution, and we have defined h k = \ * + ' w(V)dV . Thus, in order for Q(V k ) to 
converge to Q(V k ) , we need p = ( K P(V)V a ^ ^ ^ j(v a ) r v v , which becomes Eq. 



(12) by p(V) = (p mt (q,p)\ 



Appendix B. Reweighting for the density of potential-energy states 

To prove Eq. (4), we first write down the distribution density of the potential 
energy h(U) in the multicanonical ensemble: 

h(U) = ^ p S[E- U(q) - K(p)]S[U - U(q)]dqdp^j [l/Q(E)]dE . 

So 

h{U) = \Z [l [E ~ U(q)] Nf/2 l ®[E - U(q)]S[U - U(q)]dq ) [l/Cl(E)]dE 
= J + J^ S[U - U(q)]dq^j(E - U) Nf/2 ~ l @(E - U)/Cl(E)dE 
= g(U)fJ(E-U) N ' /2 - 1 /n(E)dE, 

where we have integrated out momenta in the first step, replaced U(q) by U in the 

second, and used g(U) = \ S[U - U(q)]dq in the last. Dividing both sides by the 

integral yields Eq. (4). 



Appendix C. Structure-Based protein model 



The protein model was described in ref. 31,32; we include it here for 
completeness. The model includes alpha-carbon atoms only, and builds a potential 
energy based on the geometry of a given protein structure: 

U= Yh(b-b () ) 2 +Y j k {6-e Q ) 2 + X^l-cosK^-^)]} 

angles 



bonds 



+ 



contacts 

«'<y'-3 



f \ 



r.. 
V v J 



r.. 
V v J 



dihedrals 
n=l,3 



I < 

not contacts 

i<j-3 



( V 2 

c 



r 

Vv J 



where b and b are the bond lengths of successive atoms along the chain in the given 
configuration and the reference values, respectively; similarly, 6 and O are the angles of 
three successive atoms, <p and <p are the dihedrals of four successive atoms. For non- 
bonded interactions, we distinguish pairs of contact and non-contact atoms; in the former 
case, we apply the fourth term with cr being the distance between i and j in the reference 

structure, while in the latter case, we apply the last term. We used the following 
parameters: k b = 100 , k e = 20 , = 1 , = 0.5 , and C = 4.0 A. 

We used a simpler criterion than the original one 31 ' 36 to define contact atoms: in 
the reference structure, if any two non-hydrogen atoms from two different residues have a 
distance r.. < r c , the corresponding alpha-carbon atoms are contacts. Note that Kouza et 

al. used the distances between alpha-carbon, instead of non-hydrogen, atoms to define 
contacts. The contacts produced by our criterion appeared to be more similar to those 
yielded by the CSU server . 



Figure Captions 

FIG. 1. a) The probability distributions for the energy in the microcanonical, canonical, 
and multicanonical ensembles, b) The ensemble weights in the canonical and 
microcanonical ensembles. 

FIG. 2. Multicanonical molecular dynamics on a Lennard- Jones system by the variable- 
temperature thermostats discussed in Sec. II. C. a-f) The probability distributions for the 
energy in the system for different methods. These variable-temperature thermostats 

produce flat energy distributions. Simple replacement of /3 with /3(E) in the canonical 
thermostat equations, labeled "Unmodified," does not produce flat energy distributions. 
The parameters are (a) Q i = 300 , (c) A VR = 0.01 , and (e) £ = 3 . For the Nose-Hoover 

chain variable-temperature thermostat in (a) and (b), M = 5 and Q k = 1 for k > 2 . Note 
in (b) the importance of a small time step if the thermostat mass Q i is too small. In panel 
(d), the alternative expression /3(E) = N f j(2K) E with A = 0.3 produced a non-flat 
energy distribution, characterized in Appendix A. 

FIG. 3. (Color online.) Simulations on the villin headpiece by a structure-based model, 
a) Estimated inverse temperature, /3(E) , as a function of energy for different cutoffs, r c . 
The default inverse temperature /3 d and reference inverse temperature /3 were the same, 
and they were set to 1.1 1 1 for r c = 5 A, 0.909 for r c = 6 A, 0.833 for r c = 6.5 A, 0.667 for 
r c = 1 A, and 0.5 for r c = 8 A. b) The probability distribution for the energy in the system, 
c) The values of /3(E) (solid, red), f3 E (dashed), and -logQ(is) + f3 E E (dotted, blue). 



d) The values of log g(U) (solid, red) and - log g(U) + PlU (dotted, blue). Panels (b)- 
(d) were obtained from the simulation with r c = 8 A. 

FIG. 4. (Color online.) Molecular dynamics with a flat density distribution on a 
Lennard- Jones system, a) Flat density distributions for three different variable-pressure 

barostats. b) Values of p(V) (thin lines), p* v (dashed), and G{V) = F(V) + p* r V (thick 
lines) for the three barostats. For the Langevin and MC variable-pressure barostats, 
velocity -rescaling was used as the independent thermostat with A^ VR =0.01. 
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Velocity-rescaling 
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